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ABSTRACT 

The study of photoionized environments is fundamental to many astrophysical prob- 
lems. Up to the present most photoionization codes have numerically solved the equations 
of radiative transfer by making the extreme simplifying assumption of spherical symmetry. 
Unfortunately very few real astronomical nebulae satisfy this requirement. To remedy these 
shortcomings, a self-consistent, three-dimensional radiative transfer code has been developed 
using Monte Carlo techniques. The code, Mocassin, is designed to build realistic models of 
photoionized nebulae having arbitrary geometry and density distributions, with both the stel- 
lar and diffuse radiation fields treated self-consistently. In addition, the code is capable of 
treating ones or more exciting stars located at non-central locations. 

The gaseous region is approximated by a cuboidal Cartesian grid composed of numerous 
cells. The physical conditions within each grid cell are determined by solving the thermal 
equilibrium and ionization balance equations. This requires a knowledge of the local primary 
and secondary radiation fields, which are calculated self-consistently by locally simulating 
the individual processes of ionization and recombination. The structure and the computational 
methods used in the Mocassin code are described in this paper. 

Mocassin has been benchmarked against established one-dimensional spherically sym- 
metric codes for a number of standa r d cases, as defined by the Lexin gton/Meudon photoion- 
ization workshops (Peauigno3 ll986tlFerland et alJl!995l [Peauigno t et alJl2001l) . The results 
obtained for the benchmark cases are satisfactory and are presented in this paper. A perfor- 
mance analysis has also been carried out and is discussed here. 
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1 INTRODUCTION 

Amongst the first numerical mod els for photoi o nized gaseous neb- 
ulae were t hose c alculated by iFloweJ 1 19681) . iHarringtonl 1 19681) 
and Rubin 1 1968). These early models included the basic phys- 
ical processes of ionization and recombination of hydrogen and 
helium, thermal balance and escape of the emitted photons from 
the nebula. However, the lack of reliable atomic data heavily 
limited the success of these models, as well as the fact that 
a number of important physical proces ses, such as charge ex- 
chang e and dielectronic re combination (Aldroyandi & Peauignot 
ll973t|Peauignot et alll978llStorevll98ll) . were not accounted for 
at the time. The evolution of photoionization modelling has gone 
hand in hand with advances made in atomic physics and computer 
technology. The application of photoionization models to a wider 
range of ions has been aided by the photo ionization cross-section 
calculations by Reilman & Manson] j l979T . and, more recently, the 
Opacity Project l Hummer et al. 1993). Compilations based on the 
latter's data ("e.g. IVerner & Yakovlevin"995l) . have made possible 
the inclusion of accurate photoionization cross-sections for many 



more ions in calculations. Ilvlendozal (l983) presented a compila- 
tion of radiative and collisional data for collisionally excited ul- 
traviolet, optical and infrared lines which was widely adopted, 
with some of these data still in use today, though most have been 
seperceded by more recent calculations such as the R-matrix calcu- 
lations of t he Iron Pr oject lHummeretalJll993l) and the Belfast 
group (e.g. iMcLaughlin & Bellll998t iRamsbottom et alj Il998t) . 
Currently, radiative and dielectronic recombination rates are still 
highly uncertain or unavailable for some ions; recent efforts to 
impro ve th e situat i on ha ve been reviewed by Nahar & Pradhan 
1 1999) and Nahar (2000). Most photoionization models include 
temperature-depen dent analytical fits to thes e reco mbination rates, 
such as those of lAldrovandi & Pequignoll ll973l) for radiative 
and high temperature dielectronic recombination, and those of 
iNussbaumer & Storevl i 19831) for low temperature dielectronic re- 
combination. 

Available computer power has increased enormously since 
the dawn of photoionization modelling. This has allowed more 
complex models to be built, including more ions, more frequency 
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points, more lines and more atomic levels. Nevertheless, the fun- 
damental assumption of spherical symmetry has always been re- 
tained. However, a glance at an image of any Galactic H II region 
will immediately demonstrate that these objects are neither spher- 
ically symmetric nor homogeneous. In addition, they usually have 
multiple exciting stars located at non-central positions in the neb- 
ula. By contrast, planetary nebulae (PNe) have only a single, cen- 
trally located, exciting star. However, even for PNe, spherical sym- 
metry is not a realistic assumption, as demonstrated by observa- 
tions with instruments such as the Hubble Space Telescope, which 
reveal an overwhelming variety in the shapes of planetary nebulae. 
These objects are very rarely circular in projection; a recent study 
inferred that about 50% of all known planetary nebulae are low ec- 
centricity ellipticals, while only about 10% are circular in projec- 
tion, with the remainder having more extreme elliptical or bipolar 
geometries lSokeJl997l.l200ll) . Some objects, for example the two 
young planetary nebulae He 2-47 and PN Ml -37, (also dubbed the 
starfish twins; Sahai 2000), show even more complicated geome- 
tries, with multiple lobes. Other PNe have FLIER s (fast , low ion- 
ization emitting regions: iBaTick et alJll993t Il994l Il998h . BRETS 
(bipolar, rotating, episodic jets: e.g. lLoDe^talJl993l) , a nsae, jets, 
knots , filaments, tails or multiple envelopes, (see e.g. IPerinottol 
l2000HCorradietalll999HGarc{a-SegurJl997l) . 

To our knowledge, only two three-dimensional photoioniza- 
tion codes have been develped so far, one bv lBaesgeneTai]<199rJ) 
and the other by Gruenwald, Viegas & de Broguiere (1997). The 
first code used a fixed number of equally sized cells and the on-the- 
spot approximation for the diffuse radiation field, with only the six 
more ab undant chemical elemen ts being taken into account. The 
work by iGruenwald et all ^997) improves on this by allowing a 
more flexible spatial grid and by using an iterative technique for 
the determination of the diffuse field and also by including twelve 
chemical elements in the simulations. 

Since most existing one-dimensional photoionization codes 
are based on the numerical solution of the equations of radiative 
transfer assuming spherical symmetry, their expansion to three di- 
mensions can be either very difficult or impractical, resulting in 
very large codes. The Monte Carlo approach to transfer problems 
provides a geometry-independent technique which can handle the 
radiation transport problem self-consistently. With this in mind, 
the Mocassin code (MOnte CArlo Simulations of Ionised Nebu- 
lae) was developed, in order to provide a three-dimensional mod- 
elling tool capable of dealing with asymmetric and/or inhomoge- 
neous nebulae, as well as, if required, multiple, non-centrally lo- 
cated exciting stars. 

Section [2] contains a description of the general Mocassin ar- 
chitecture and of some of the main computational methods used 
in the code. The code has been benchmarked against established 
spherically symmetric one-dimensional photoionization codes for 
a set of standard nebulae and in Section [3] we present the results 
of this benchmarking, together with a performance analysis of the 
codes. In section|4]we discuss the results of the benchmarking and 
present some general guidelines on how to run the code efficiently. 



2 DESCRIPTION OF THE MONTE CARLO CODE 
2.1 Background 

The Monte Carlo method has been widely applied to a variety of 
astrophysical problems, such as the penetration of ultraviolet ra- 
diation into the interiors of uniform or lumpy, interstellar clouds 



iFlannerv et alll98(j;lBoisse1l990l). resonance-like scattering in ac- 
cretion disc winds <KniggeetalJll995l) and polarizatio n maps for 
the circumstellar envelopes of protostars iFischer et alJll994h . In 
the examples described above the absorption and scattering coef- 
ficients are not coupled to the radiation field and, therefore, these 
problems do not require solution by iteration. 

However, Monte Carlo techniques have also been used for 
dust radiative equili b rium calculati o ns for some time, see e.g. 
Lefevre et al] ll982l) . Lefevre et all 1983) and, more recently, 
Wolf et al. 1 1999). These authors use a technique in which stellar 
and diffuse photon packets are emitted separately; the number of 
diffuse photon packets (i.e. packets emitted by the dust) is deter- 
mined by the dust grain temperature, which in turn is determined 
by the balance between the number of absorbed and emitted photon 
packets. An initial guess for the dust grain temperature is provided 
by the number of packets absorbed, and the iteration continues un- 
til the grain temperatures converge. Using this method the stellar 
luminosity is not automatically conserved during the Monte Carlo 
simulation; only after the grain temperatures have reached conver- 
gence is the stellar luminosity approximately conserved. The con- 
vergence of such codes is often very slow and requires a large num- 
ber of iterations and simulation quanta in order to reach the required 
accu racy. 

iBiorkman & Woocl 1200 ll) have described a general radiative 
equilibrium and temperature correction procedure for use in Monte 
Carlo radiative transfer codes having sources of temperature- 
independent opacity, such as dust. Their technique makes use of 
information naturally given by the Monte Carlo method, which, by 
tracking every photon/energy packet, makes it easy to determine 
where in the simulation grid energy is being absorbed. When en- 
ergy is deposited at a given location, following a packet's absorp- 
tion, the local medium is heated. Whenever this occurs the new 
local temperataure is calculated and the packet is then re-emitted 
accordingly. The packets are followed in their path through the re- 
gion, as they undergo scatterings and absorptions followed by re- 
emissions, with the temperatures being updated after each event, 
until the packets reach the edge of the nebula and escape to infinity, 
hence contributing to the emergent spectrum. Once all the stellar 
photon packets have escaped, the resulting envelope temperature 
and the emergent spectrum are correct without the need of any fur- 
ther iterations. 

A great limitation of Bjorkman & Wood's method is that it 
cannot be applied to situations where the opacities are temperature- 
dependent, as is the case in photoionized nebulae. There are two 
reasons for the failure of this method when the opacity varies with 
the local temperature: firstly, the number of photon packets ab- 
sorbed by the cell prior to a temperature update would be either too 
small or too large, and, secondly, a change in temperature would 
also imply a change of the interaction locations of previous pack- 
ets, signifying that the paths of the previous photon packets should 
have been different. While, it is clear that, when dealing with pho- 
toionised gas, Bjorkman & Wood's technique is not applicable, 
their work is nevertheless very enlightening and should be taken 
into account for further developments of the Mocassin code, when 
a treatment for dust grains will be introduced. 

A recent example of the application of the Monte Carlo tech- 
nique to pro blems requiring solution by iteration is the work of 
lLucvliT999T) . who obtained the temperature stratification and emer- 
gent spectrum of a non-grey spherically symmetric extended stel- 
lar atmosphere in LTE. H i s resu lts show very good agreement with 
the predictions of lCastoJ ll974l) . hence demonstrating the validity 
of the Monte Carlo techniques applied, some of which were also 



Mocassin 3 



used in the development of Mocassin. The current work folows the 
approach described by iLucyj < 19991) and also applied in the one- 
dimensional code developed by Och, Lucy & Rosa (1998). They 
employed a different Monte Carlo treatment of the radiative trans- 
fer in order to iteratively determine the temperature and ionisation 
stratification for a spherically symmetric photoionised nebula of 
uniform density. So me of the techniques that they used are also 
described in detail bv lLucvl 1 1999, 20011 120021) . The basic concept 
is that, when calculating radiative equilibrium temperatures, con- 
servation of stellar luminosity is more important than the details 
of the spectral energy distribution. With this in mind conservation 
of stellar luminosity is enforced by using energy packets of con- 
stant net energy throughout the simulations. Moreover, all absorbed 
packets are re-emitted immediately after every absorption event. 
The frequencies of the re-emitted energy packets are determined 
by the local gas emissivities. Although the frequency distribution 
of the re-emitted packets will not be correct until the nebular tem- 
peratures have converged, this method naturally enforces radiative 
equlibrium at each point in the nebula and so naturally provides 
conservation of energy. This not only results in a simpler co de but 
also makes the convergence of the gas temperatures easier iLucvl 
1 19991 2001). Energy packets will be discussed in more detail in 
section 12.21 



2.2 Energy Packets 

The main principle of our treatment of a photoionized nebula con- 
sists of locally simulating the individual processes of ionization and 
recombination. The radiation field is therefore expressed in terms 
of energy packets, e(y), which are the calculation quanta. e(v) is a 
packet consisting of n photons of frequency v such that 



e[y) — nhv 



(1) 



In addition, we take all packets to have constant energy eo. There 
are several reasons for choosing to work with monochromatic, in- 
divisible packets of radiant energy instead of photons. First of all, 
energy packets are more computationally economic and, also, since 
they all have the same energy, then those packets emitted in the in- 
frared will contain a larger number of photons wh ich, as a conse- 
quenc e, will not have to be followed individually lAbbott & Lucvl 
1985). Note that all energy packets are followed until they escape 
the nebula, including infrared energy packets. This is in order to 
allow the introduction of dust particles into the radiative transfer 
treatment of Mocassin, which is planned for the near future. Also, 
as the total stellar luminosity, L t , is evenly split amongst the stellar 
energy packets, the energy carried by a single packet in the time 
interval A t, which represents the duration of the Monte Carlo ex- 
periment, is given by 



~N 



eo 
At 



(2) 



where N is the n umber of energy packets used in the simulation 
lQchetal] fl998). Most importantly, the use of constant energy 
packets is a natural way of imposing strict energy conservation at 
any point in the nebula iLucvll999T) . So, when a packet of radiant 
energy e(u a ) = eo is absorbed, it is immediately re-emitted with a 
frequency v c , which is determined according to a frequency distri- 
bution set by the gas emissivity of the current cell. The packet emit- 
ted, e(v c ), will then have the same energy as the absorbed packet, 
e(v a ), meaning that only the number, n, of photons contained in 
the packet is changed. 



2.3 Initiation 

In our modelling the gaseous region is approximated by a three- 
dimensional Cartesian grid, where the ionising source can be 
placed at the centre of the grid or anywhere else in the grid. This 
feature is very useful when dealing with axisymmetric nebulae, 
since, by placing the source in a corner of the grid, we need only 
consider one eighth of the nebula, which can then be reconstructed 
in full at the end of the simulation. This allows the running of mod- 
els with much higher spatial resolution than those which would be 
possible if a full nebula had to be considered, by putting the source 
in the centre and, therefore, not making use of any symmetry prop- 
erties of the object. Switches built inside the code allow the user 
to specify whether the nebula has some degree of symmetry and, if 
so, whether the symmetry is to be used. 

Inside each grid cell all nebular properties, such as the mass 
density of the gas, p; the electron temperature and density, T c and 
N e ; and the frequency dependent gas opacity and emissivity, k„ 
and j u , are constant by definition. Thermal balance and ionisation 
equlibrium are imposed in each grid cell in order to obtain the phys- 
ical conditions in the local gas. 

The energy packets are created at the position of the ionis- 
ing source and they all carry the same energy eo, as discussed in 
the previous section. The frequency, v, of each individual packet 
emitted is derived from the input spectrum of the ionising source 
according to the probability density function 

F v dv F v dv 



p{v) 



(3) 



where F v is the stellar flux and is the stellar radius. This is then 
the probability of an energy packet being emitted with a frequency 
lying in the interval (u, v + dp). The upper and lower integration 
limits, v m in and u max , have to be chosen properly, depending on 
the input spectrum, in order to ensure that the bulk of the radia- 
tion is included in the frequency range. As the source emits energy 
isotropically, the direction of travel of every energy packet emitted 
is chosen randomly. This is done by choosing two random num- 
bers, a and /3, in the interval [0, 1], and calculating the following 
quantities: 

w — 2a — 1 



t = \J\ - w 2 

6 = 7r(2/3-l) 

u = tcos8 

v — tsinO 



(4) 



The rando m unit vector in Cartesian coordinates is then (u, v, w) 
jHarries & Howarthl 199% . 



2.4 Trajectories 

Once a stellar packet is created at the source and launched into 
the nebula, its trajectory must be computed as it undergoes absorp- 
tions followed by re-emissions due to bound-free and free-free pro- 
cesses. The trajectory ends when the packet reaches the edge of the 
nebula, where it escapes to infinity and contributes to the emergent 
spectrum. 

There are two methods to track the packets and determine the 
locations of the absorption events. Consider a packet of frequency 
v p , emitted in the direction u. The first of these methods consists 
of calculating the run of optical depth, t„ p , at the energy packets' 
frequency v v , from the location at which the packet is emitted to 



4 Ercolano et al. 



the edge of the ionised region along the direction of travel, u. The 
probability of absorption along that path is then given by 

p(7V„) = e^p (5) 

and the normalised cumulative probability function is given by 



P(l) = 



(6) 



where t v (V) is the optical depth to the absorption event and I is 
the path length. The position at which the energy packet will be 
absorbed will then be determined by choosing a random number 
in the interval [0, 1] and comparing it against P(l). In reality, it 
is more convenient to use the inverse approach, where the optical 
depth from the energy packet source to the event can be derived 
from the inverse of equation [5] 



(7) 



where Ur is a random number in the interval [0, 1]. Once r v (I) 
has been ca lculated then the path length can be directly derived 

lHarries & Howarthl 19971) . 

The second method was sugg estedbv lLucvl fl999) and it con- 
sists of testing whether an absorption event occurs, on a cell by 
cell basis. In other words, assume that, within each uniform cell, 
the random path of a packet between events is given by equation^ 
which corresponds to a physical displacement, I, given by 



r Vp = K„pl 



(8) 



where k„ and p are the frequency dependent absorption coefficients 
and the density of the current cell respectively. The method then 
consists of checking whether the displacement I is large enough 
to carry the packet out of its current cell. If this is the case, the 
packet is moved along its direction of travel, u, up to the boundary 
with the adjacent cell, where a new value for Ur is cast, giving a 
new t v , and any further movement of the packet in this new cell 
is to be followed. Alternatively, if the displacement I is not large 
enough to carry the energy packet across the next boundary, the 
packet will be absorbed and then re-emitted at the end-point of the 
displacement. Lucy also clarifies in his paper that the selection of a 
new value of t„ p at the crossing of a boundary does not introduce 
a bias since a photon always has an expected path length to its next 
event corresponding to t v = 1, regardless of the distance it might 
already have travelled. 

In this work both methods were implemented in the code, in 
turn, in order to test their respective performances. The first method 
proved to be much more computationally expensive then the sec- 
ond. This is due to the fact that, in order to track down the posi- 
tion at which an energy packet is absorbed, using our knowledge 
of T Vp (l), an array searching routine has to be used to locate the 
index of t„ (i) within the array of optical depths calculated from 
the packet's source to the edge of the nebula. Although the search- 
ing procedure employs a bisection technique, which makes it quite 
efficient, the large number of calls to it, due to the large number of 
energy packet interactions within a simulation, means that nearly 
60% of the run time is spent carrying out these searches. The sec- 
ond method does not require any calls to the array searching rou- 
tine, as the packets are followed step by step through the nebula, 
and this results in the run time being considerably reduced. The 
current version of Mocassin therefore uses Lucy's appoach to track 
the energy packets throughout the nebula. 

Finally, the direction of travel of the newly emitted diffuse 



energy packets (i.e. those packets re-emitted immediately after an 
absorption event) needs to be determined. Since absorption and re- 
emission are two independent events, the diffuse packets are emit- 
ted isotropically and therefore their direction of travel is chosen 
randomly using equations|4] 



2.5 The Mean Intensity 

The success of a Monte Carlo model often relies on the careful 
choice of appropriate estimators. Monte Carlo estimators provide 
the means to relate the quantities we observe during our Monte 
Carlo experiment to the physical quantities we want to determine. 
In a photoionization model, a measure of the radiation field is 
needed, namely th e mean intensity, J v . 

In the work of O ch et alJll998l) . the Monte Carlo estimator of 
J v is constructed by using the definition of the specific intensity, 
/„, in spherical coordinates, (r, 8), as a starting point: 



AE = I v (r,9)AA | cos6 | A vAujAt 



(9) 



where A A is the reference surface element, 8 is the angle between 
the direction of light propagation and the normal to the surface A A 
and A u> is the solid angle. The mean intensity can then be obtained 
from this by calculating the zero order moment of /„, which gives 



4nJ v (r) = / l v doj = — — 



ae^A 



i 



i 



A t l—' cosft A A A v 



(10) 



by comparison with equation|5| The sum is over all packets iVk with 
frequency lying in the interval (y, v+d v), crossing A A at an angle 
9. As discussed above, A E/A t represents the energy carried by a 
single packet in the time interval A t, since AE — eo, which is 
given by equation|2| Equation llOl then provides a relation between 
the Monte Carlo observables (i.e. the number of energy packets 
with frequency lying in the interval (y, v + dv), crossing A A at 
angle 8 and the mean intensity of the radiation field, which is the 
required physical quantity. 

The use of Och et al.'s estimators for J„, however, becomes 
problematic in the non-spherically symmetric case, since the refer- 
ence surface for the volume elements in an arbitrary two- or three- 
dimensional coordinate system might not be unique or as obvious 
as in the one-dimensional case. In our work, a more general ex- 
pression for the es timator of J v is sought, and, therefore, following 
Lucy's argument <Lucvlll999l) . an estimator for J v is constructed 
starting from the result that the energy density of the radiation field 
in the frequency interval {y, u + dv) is A-Kj v du/c. At any given 
time, a packet contributes energy e(u) = eo to the volume element 
which contains it. Let I be a packet's path length between succes- 
sive events, where the crossing of cell boundaries is also considered 
an event; the contribution to the time averaged energy content of a 
volume element, due to the I fragments of trajectory, is eoSt/ At, 
where St — l/c. From this argument it follows that the estimator 
for the volume element's energy density can be written as 



AnJ v dv eo 1 



I 



c At V ^ c 

d v 



(ID 



where V is the volume of the current grid cell and the summation 
is over all the fragments of trajectory, I, in V, for packets with fre- 
quencies lying in the interval (y, v+d v). Again, a relation between 
Monte Carlo observables (i.e. the flight segments, /) and the mean 
intensity of the radiation field, J„ has been obtained. Moreover, 
equation \H] is completely independent of the coordinate system 
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used and, indeed, of the shapes of the volume elements, V . An- 
other important aspect of this approach is that all packets passing 
through a given grid cell contribute to the local radiation field even 
without being absorbed; this means that equation fTTI returns estima- 
tors of the radiation field even in the extremely optically thin case 
when all packets pass through the nebula without any absorption 
events. From this argument it follows that this technique allows a 
much better sampling and, hence, in general, much less noisy re- 
sults compared to other techniques based on estimators for which 
only packets absorbed within a given volume element count. 

2.6 Gas emissivity and the diffusion of energy packets 

As we have already discussed in previous sections, after an energy 
packet is absorbed, a new packet is re-emitted from the same loca- 
tion in a random direction. The frequency of the re-emitted packet 
is calculated by sampling the spectral distribution of the total lo- 
cal emissivity, j* oi . In order to satisfy the thermal balance implied 
by the Monte Carlo model, all major emission processes have to 
be taken into account, including the complete non-ionizing nebu- 
lar continuum and line emission, since they are part of the energy 
budget. The non-ionizing radiation generated in the nebula is as- 
sumed to escape without further interaction and constitutes the ob- 
servable spectrum which can then be compared with observations. 
The following paragraphs are concerned with the description of the 
individual contributions to the total emissivity. 

The continuum emission due to H I, He I, He II and to heavier 
ions is included. The H I continuum can be divided into the Ly- 
man continuum, which is capable of ionizing H, and the Balmer, 
Paschen, etc. continua, which are not capable of ionizing H. The 
emissivity in the Lyman continuum is calculated directly from a 
combination of the Saha and Milne relations: 



h 2 



\3/2 



C OJi+1 



a„(X l 



-h(v-v )/kT m Y i+1 



X I+1 Ar o (12) 



where u>i and Wi+i are the ground state statistical weight of the ions 
involved, X l+1 is the abundance of the recombining ion, a v (X.') 
is the photoionization cross section and vo is the photoionization 
threshold. The emissivity of the other series continua are obtained 
by interpolation of published data (Ferland 198C). A similar ap- 
proach is used for the He I and the He II continua, where for 
frequencies greater than 1.8 Ryd and 4.0 Ryd, respectively, equa- 
tion^] is used, and the emissivities at lower freque ncies are ob- 
tained by interpolation of the data published by Brown & Matthews 
1 1970) for the He I series and bv lFerland ll98(J) for the He II se- 
ries. The continuum emissivity of heavy elements is also calculated 
using equation fL?! In the hydrogenic case (i.e. H I and He II), the 
tw o-photon continuum is calculated using the f ormalism described 
bv lNussbaumer & SchmutJ Il984h : the data of iDrake et"aT] 1 19691) 
are used for He I. Recombination lines between lower levels n=2 
through 8 and upper levels n=3 through 15, for H I, and lower lev- 
els n=2 through 16 and upper levels n=3 through 30 for He II, are 
calculated as a function of temperature according to the case B data 
published by IStorev & Hummed ll995l) . The He I recombination 
lines are calculated as a function of temperature using the data of 
lBeniaminetai]jl999l) . In general, He I singlet lines follow Case B 
whereas triplet lines follow Case A (as there is no n = 1 level for the 
triplets). Transitions to the 1 X S ground state of He I produce lines 
which are capable of ionizing H and low ionization stages of higher 
elements. In particular, the emissivities of the He I Lyman lines 
from n=2 through n=5 lBrocklehurstlll972l) and the intercombina- 
tion lines corresponding to the transitions 2 3 S-1 1 S and 2 3 P-1 1 S are 



estima ted as a function of temperature using the data of [Robbins 
1 1968). The contributions due to these lines to the total energy dis- 
tribution, from which the probability density functions are derived, 
are added into the respective energy bins. Similarly, He II Lyman 
lines can ionize both neutral hydrogen and neutral helium, as well 
as some of the low ions of heavier elements. Therefore the emissiv- 
ities of He II Lyman lines with upper levels from n=2 through n=5 
(fits to Storev & Hummer 1995) are also estimated as a function 
of temperature and their contributions to the total energy distribu- 
tion added into the respective frequency bin, as for the He I lines. 
This method is based on the fact that all emission profiles are cur- 
rently treated as S functions and the line opacity is assumed to be 
zero; and the absorption of energy packets is only due to the contin- 
uum opacity. Finally, the emissivities of the collisional lines of the 
heavier ions are calculated. This is done by using matrix inversion 
procedures in order to calculate the level populations of the ions. 
Appendix 1 contains references for the atomic data used for each 
ion. 

The energy distribution is derived from the total emissivity, 
summing over all the contributions in a particular frequency in- 
terval. The non-ionizing line emission is treated separately, since, 
whenever such line packets are created, they escape without further 
interaction . 

Once the line and continuum emissivities have been calcu- 
lated, the probability that the absorption of an ionizing energy 
packet will be followed by the emission of a non-ionizing packet 
is given by: 



p - 

J esc — 



(13) 



where f m „ is the higher limit of the frequency grid; the j xi are the 
emissivities of the non-ionizing recombination lines of all species 
considered; j„ is the frequency dependent continuum emissivity; 
j^ei an d JHeii are me contributions due to those recombination 
lines of He I and He II which are capable of ionizing neutral hy- 
drogen and neutral helium. The choice between the re-emission of 
an ionizing photon or a non-ionizing one is made at this point in 
the code. 

If an ionizing energy packet is to be re-emitted, then the new 
frequency will be calculated according to the normalised cumula- 
tive probability density function for the ionizing radiation, given 

by 



p{v) = 



Q^+EiHei^EiU 



(14) 



where, as usual, the contributions due to the He I and He II lines 
are added in the corresponding frequency bins. If a non-ionizing 
energy packet is to be re-emitted, then its frequency must be deter- 
mined from the probability density function for non-ionizing radia- 
tive energy, which is analogous to eauation fT4l 

2.7 The Iterative Procedure 

An initial guess of the physical conditions in the nebular cells, such 
as the ionization structure, electron temperature and electron den- 

1 Resonance lines longward of 912 A (e.g. C IVAA1548, 1550) may, in fact, 
diffuse out of the nebula via resonant scattering and may also be absorbed 
by dust during such diffusion. A treatment of dust grains will be included 
in future developments of the Mocassin code, and such effects may then be 
accounted for. 
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sity, needs to be specified before the simulation can begin. Proce- 
dures in Mocassin have been constructed such that only an initial 
guess at the electron temperature (which is initially set to a constant 
value throughout the nebula) must be included in the input file. Mo- 
cassin can then guess an initial ionization structure and, hence, the 
electron density. However, if the output of a one dimensional model 
(or a combination of more than one of them) is available, there are 
also procedures built into Mocassin to map these onto the three- 
dimensional Cartesian grid, by using simple interpolation routines. 
A one-dimensional mode option was implemented in Mocassin for 
this purpose. Several tests have shown that while the choice of the 
initial conditions has, of course, no influence on the final result of 
the simulation, it can, however, have an inpact on the number of 
iterations required to reach convergence. It is hard to quantify the 
number of iterations required for convergence by each method, in 
particular, it depends strongly on the initial temperature input used 
in the first method, and, when applying the second method, on the 
deviation of the actual three-dimensional geometry from the sim- 
plified one-dimensional model used. However, with sufficient en- 
ergy packets, the benchmark models described here should be fully 
converged in approximately fifteen to twenty iterations. A strategy 
to speed up the simulations is described in Section 3.1. 

Once the initial conditions are specified, the frequency depen- 
dent total emissivities are calculated in each grid cell in order to 
set up the probability density functions for re-emitted radiation, 
which are used for the determination of the frequency distribution 
of the re-emitted energy-packets during the Monte Carlo simula- 
tion. The energy packets are then fired through the grid and their 
trajectories computed. Once all the energy packet trajectories have 
been computed, the Monte Carlo estimators for the mean inten- 
sity of the stellar and the diffuse radiation fields can be obtained, 
as described in Section [23] The ionization fraction and the elec- 
tron temperatures and densities must now be updated to be self- 
consistent with the current estimates of the radiation field at each 
grid point. This means solving the local ionization balance and ther- 
mal equilibrium equations simultaneously. The entire procedure is 
repeated until convergenge is achieved. The convergence criterion 
that is used in this work is based on the change of the local hy- 
drogen ionization structure between successive iterations. In some 
cases, however, this is not a suitable convergence criterion (e.g. in 
hydrogen-deficient environments), for this reason, other criteria are 
also implemented in the code (e.g. based on the change of the lo- 
cal helium ionization structure, or of the local electron temperature 
between successive iterations), and these can be easily selected by 
using the appropriate switches in the input file. 

2.8 Comparison of the Model with Observations 

When the model has converged to its final solution, the output spec- 
trum can be computed and compared with the results obtained from 
other models or with observational data. The total luminosity of 
the nebula emitted in various emission lines longward of the Ly- 
man limit can be obtained by using two methods. The first method, 
which is only available to Monte Carlo codes, consists of summing 
up the number of energy packets in the given line, Nu ne , over the 
grid cells. Hence, the power emitted in the line is given by 

imax jmax k max 

Lline = -^T ^ ^ N line{Xi,Vj, Zfc) (15) 

i=l j=l fc=l 

where ^ is given by equation [2] The second method consists of 
using the values of the local electron temperature and ionic abun- 



dances given by the final converged model solution to obtain the 
line emissivities for each grid cell. The luminosity of the nebula in 
any given line can then be calculated easily by summing the emis- 
sivity of the required line over the volume of the nebula. 

A comparison of the results obtained using the two methods 
described above, provides an indication of the level of accuracy 
achieved during the simulation, as the two methods will give con- 
sistent results only if enough energy packets have been used to yield 
good statistics for every line. In general, the second method (for- 
mal solution) yields the most accurate results, particularly for weak 
lines, which may emit relatively few photons. For the benchmark 
cases presented here, reasonable accuracy was deemed to have been 
achieved when the fluxes of the strongest transitions obtained using 
the pure Monte Carlo method were within 10% of those obtained 
using the formal solution. Both methods can also be used to cal- 
culate line of sight results and to simulate long-slit observations. 
However, just as for the calculation of the integrated spectrum, the 
formal solution method is to be preferred, as it yields the most ac- 
curate results, particularly for the weaker lines. 

In addition to the integrated emergent spectrum, other useful 
comparisons with the observations can be carried out, e.g. projected 
images of the final model nebula in a given line or at a given con- 
tinuum frequency can be produced for arbitrary viewing angles. 
These can be compared directly with nebular images obtained in 
an appropriate filter. Mocassin computes and stores the physical 
properties of the nebula, as well as the emissivities of the gas at 
each grid point; these can be fed into IDL plotting routines in or- 
der to produce maps (Morisset et al., 2000). Also, by assuming a 
velocity field, line spectral profiles can be produced, together with 
position-velocity diagrams. These can be compared with observa- 
tions, if available, to deduce spatio-kinematic information about the 
object being studied. More information about the original IDL rou- 
tines is given bv lMorisset et alj i20obT) and lMonteiro et alJ (2000). 
Details of the actual application to Mocassin's grid files are avail- 
able in a compan ion paper on the modelling of the planetary nebula 
NGC 3918 lErcolano et alfeOOl Paper II). 

At the end of each Monte Carlo iteration the physical quan- 
tities which characterise the grid are written out to disc into three 
files, namely gridl.out, grid2.out and grid3.out. The first file con- 
tains the local electron temperature and density as well as the gas 
density at each grid cell, the second the ionization structure of the 
nebula and the third a number of model parameters, including the 
number of energy packet to be used in the simulation. These files 
are used in conjuction with a warm start driver, which allows an in- 
terrupted simulation to be resumed from the end of the last Monte 
Carlo simulation. This means that once a simulation has been inter- 
rupted the number of energy packets used (and indeed other model 
parameters, if required) can be adjusted, before the simulation is 
restarted, by modifying the file grid3.out. This feature can be used 
to speed up the simulations by using the following approach. The 
first few iterations are run using a lower number of energy packets 
than actually needed; so, for example, if the optimum number of en- 
ergy packets for a given model is 10 6 , then the first few iterations 
can be carried out using only 10 5 packets, hence reducing the run 
time for these by a factor of ten. This will result in about 50%-60% 
of the grid cells converging; in general, the inner cells converge 
more quickly, due to the larger number of sampling units available 
there (due mainly to geometrical dilution and to the reprocessing 
of energy packets to non-ionizing energy packets). At this point the 
simulation is interrupted and then resumed, after having adjusted 
the number of energy packets to the final required value (i.e. 10 6 , 
in the previous example). Final convergence will then be achieved, 
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in most cases, within four or five further iterations. The actual num- 
ber of iterations required depends on the number of energy packets 
used: the larger the number of sampling quanta available at each 
cell, the quicker the cells will converge to a solution. The numbers 
quoted above, however, also depend on each particular model's ge- 
ometry and optical thickness. 

2.9 General Architecture 

The Mocassin code was written using the Fortran 90 programming 
language. The code was developed and run initially on a Com- 
paq(Dec) XP1000 with a 500 MHz CPU and 1 Gb of memory and 
a preliminary serial version of the code still exists. A fully paral- 
lel version of the code has since been developed using Multiple 
Processes Interface (MPI) routines and it currently runs on a Sil- 
icon Graphics Origin 2000 machine with 24 processors and 6 Gb 
of memory and a SUN Microsystems Sunfire V880 machine with 
16 processors and 64 Gb of memory. Monte Carlo simulations are, 
by their nature, very parallelizable problems and, indeed, Mocassin 
can achieve a linear speed-up, i.e. a speed-up that is directly pro- 
portional to the number of processors used. A detailed description 
of all the Mocassin modules, input commands and output files is 
given bv lErcolandfeOOl PhD Thesis). A copy of the code is avail- 
able from the author (be@star.ucl.ac.uk) together with the relevant 
thesis chapters. 



3 APPLICATION TO BENCHMARK CASES 

Numerical simulations of photoionized nebulae are very complex 
and a number of factors, such as numerical approximations and as- 
sumptions, and the complexity of the calculation itself, introduce 
a degree of uncertainty into the results. For this reason, it is im- 
portant for modelers to have certain standards of comparison, in 
order to identify problems in their codes and to reach an adequate 
degree of accuracy in their calculation. A series of meetin g have 
been held, beginning in Meudon, France, in 1985 I Peauienot 1986) 
and in Lexington, Kentucky , first in 19 95 iFerland et alll995T) and 
again in 2000 (Peauienot et al. 2001), in order to define a set of 
benchmark cases which could be used by all photoionization mod- 
elers to test their codes against. The benchmarks which resulted 
from these meetings include H II regions, planetary nebulae, nar- 
row line regions (NLRs) of AGNs and X-ray slabs. Mocassin does 
not have, at present, the capability to treat NLRs and X-ray slabs, 
as some relevant physical processes, such as Compton heating and 
inner shell ionization, are not yet included. For this reason, only the 
H II region and planetary nebula benchmarks are performed in this 
work. The expansion of the code to include high energy processes 
is planned in the future. 

Results from several other codes are available for comparison; 
these are all one-dimensional codes and, apart from differences in 
the atomic data used by each of them, their main differences lie 
in the treatment of the diffuse radiation fiel d transfer. A b rief de- 
scription of each of these codes is given bv lFerland et alj ll995l) . 
Although the majority of these codes have evolved somewhat since 
the 1995 Lexington meeting, mostly via the updating of the atomic 
data sets and the inclusion of more and specialised physical pro- 
cesses, their basic structures have stayed the same. The seven codes 
included for comparison are G. Ferland's Cloudy (GF), J.P Harring- 
ton's code (PH), D. Pequignot's Nebu (DP), T. Kallman's XStar 
(TK), H. Netzer's Ion (HN), R. Sutherland's Mappings (RS) and 
R. Rubin's Nebula (RR). Only two of these codes, the Harrington 



Table 1. Lexington 2000 benchmark model input parameters. 



Parameter 


HII40 


HII20 


PN150 


PN75 


L(BB)/10 37 (f^§) 


308.2 


600.5 


3.607 


1.913 


T(BB)/10 3 K 


40 


20 


150 


75 


R in /10 17 cm 


30 


30 


1 


1.5 


mr/cm - 3 


100 


100 


3000 


500 


He/H 


0.10 


0.10 


0.10 


0.10 


C/Hx 10 5 


22. 


22. 


30. 


20. 


N/Hx I0 5 


4. 


4. 


10. 


6. 


O/Hx 10 5 


33. 


33. 


60. 


30. 


Ne/Hx 10 5 


5. 


5. 


15. 


6. 


Mg/Hx 10 5 






3. 


1. 


Si/Hx I0 5 






3. 


1. 


S/Hx 10 5 


0.9 


0.9 


1.5 


1. 



Elemental abundances are by number with respect to H. 



code and Rubin's Nebula, treat the diffuse radiative transfer exactly. 
The others use some versions of the outward-only approximation of 
varying sophistication. In this approximation all diffuse radiation is 
assumed to be emitted isotropically into the outward half of space. 

The predicted line fluxes from each code for each benchmark 
case are listed in Tables|4]to[7] together with the volume-averaged 
mean electron temperature, weighted by the proton and electron 
densities, N p N e , <T[N P ,N C ] >, the electron temperature at the 
inner edge of the nebula, Tinner, and the mean ratio of fractional 
He + to fractional H + , -r , which represents the fraction of 
helium in the H + region that is singly ionized. <T[N H+ ,N e ] > 



and 



<He+> 



are ca lculated according to the following equations 



<T[N P ,N C ] >-- 



j N C N P T C dV 
J N e N p dV 



(16) 



and 

< He+ > 



(17) 



w(H) / N c N Ha+ dV 

< H+ > n(He) J N e N p dV 

where N c and N p are the local electron and proton densities, 
respectively, N Hc + is the density of He + , and n(H) and 7i(He) are 
the total hydrogen and helium densities. 

TableQlists the input parameters for all the benchmark models 
dicussed here. All the benchmark cases listed in TableQwere cal- 
culated using both the three-dimensional and the one-dimensional 
mode of Mocassin and both sets of results are included here for 
comparison. It is clear from Tables [4] to Q that the results of the 
three-dimensional and one-dimensional modes of Mocassin are 
consistent with each other. The small differences that do exist can 
be entirely attributed to the coarseness of the grids used for the 
three-dimensional calculations. The aim of the benchmarking de- 
scribed in this work is to assess the reliability of Mocassin in its 
fully three-dimensional mode, for this reason the one-dimensional 
mode results will not be included in the following performance 
analysis; moreover the inclusion of two sets of results from what 
is, essentially, the same code would introduce a bias in the median 
and isolation factors calculations described below. Finally, to avoid 
any confusion, any mention of Mocassin throughout the rest of this 
paper refers to the fully three-dimensional version of the code, un- 
less otherwisestated. 

Figures Q and |2| show the electron temperatures (top panels) 
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and the fractional ionic abundances of oxygen (middle panels) and 
carbon (bottom panels) for the four benchmark cases analysed. The 
ionic abundances in every cell in the ionized region are plotted 
against radial distance from the star. These plots are interesting not 
only because they provide a clear picture of the overall temperature 
and ionization structure of each model, but also because from the 
scatter of the data points one can estimate the accuracy of the final 
results. (Note that such plots are only meaningful in the spherically 
symmetric case.) 

Four benchmark model nebulae were computed, two H II re- 
gions and two planetary nebulae. These benchmarks were designed 
to be uncomplicated yet to test different aspects of the modelling 
(see lFerland et alll995h . The nebulae are homogeneous in density 
and, for simplicity, blackbodies are used as the ionizing sources 
instead of model stellar atmopheres. 

F ollowing the analysis of Pequignot (see IPeauignot et alJ 
l200ll) . isolation factors, if's, were computed for each predicted 
quantity in each case study. These are defined as the ratio of the 
largest to the penultimate largest value of a given output quantity 
or the ratio of the penultimate smallest value to the smallest value. 
These ratios are computed with the intention to identify aberrant 
values. A large if can be attributed to a number of factors, but of- 
ten these can be attributed to a difference in the atomic data used 
by each modeler. A list of the number of if's larger than 1.01, 
1.03, 1.10, 1.30 and 2.00 is given in Table|8| for each benchmark. 
After analysing the benchmark results obtained b y all the m odel- 
ers wh o participated in the Lexington workshop. IPeauignot et alJ 
1200 ll) suggested that an isolation factor larger than 1.30 is indica- 
tive of a significant departure and a possible problem. A large num- 
ber of occurences of if's > 1.30 should either have an acceptable 
explanation or lead to corrections to the code. 

The number of results not predict ed by any given code is listed 
in the No pred row of Table|S| Pequign ot et al Ji200lh also noted, in 
the proceedings of the November 2000 Lexington meeting, that the 
lack of a prediction for a particular observable may simply reflect a 
lack of interest by the modeller in it; on the other hand, a frequent 
occurence of No pred may also indicate limitations in the predictive 
power of a given c ode. 

As argued bv lPeauignot et all 1200 ll) . a large error can be in- 
troduced when the average over a small sample containing a num- 
ber of aberrant values is taken. In order to minimise this error, me- 
dian values are calculated instead of averages and these are given 
for each observable listed in Tables|4]to0 in the column labelled 
Med. The medians are calculated to the precision shown in Tables|4] 
toQ Table|5|lists the number of median values scored by each code 
for each benchmark, i.e. the number of times the code was the clos- 
est to the median value. When a median value is shared by two or 
more codes the score is given to each one, therefore the sum of the 
median values scored by all the codes is higher than the number of 
observables (the column labelled Total in Table[9). 

3.1 Sampling Requirements 

Table|2|lists the optical depths at the ionization threshold frequen- 
cies for H° , He and He + , at the outer edge of the grids, for the four 
benchmark models analyzed here. For each model, the number of 
grid points is also given (column 5), together with the number of 
energy packets used, -/Vpackets, according to the two-step stategy 
described above, first to achieve convergence in 50%-60% of the 
total number of grid cells (> 50%, column 6) and then to achieve 
total convergence (> 95%, column 7). Table|2|shows that the softer 
the ionizing radiation field, the larger the number of energy packets 



Table 2. Summary of the number of energy packets needed for > 50% and 
> 95% convergence (see text for explanation) for each of the benchmark 
cases 



Case 


H° 


7~edge 






-Vpackets 




He 


He+ 




>50% 


>95% 


HII40 


4.79 


1.15 


177.8 


13x13x13 


5-10 5 


5-10 6 


HII20 


2.95 


1.13 


91.2 


13x13x13 


5-10 6 


5-10 7 


PN150 


34.0 


6.87 


57.9 


13x13x13 


310 s 


3-10 6 


PN75 


1.16 


0.24 


31.5 


13x13x13 


4-10 5 


4-10 6 



Table 3. Deviation of the Monte Carlo method from the formal solution for 
the prediction of some significant line fluxes in the benchmark models. 



Line 


HII40 


HII20 


PN150 


PN75 


H/3 


2.7% 


9.5% 


5.8% 


2.8% 


He I 5876 A 


5.2% 


6.3% 


0.96% 


4.5% 


[N II] 6584 A 


7.6% 


4.9% 


8.5% 


4.8% 


[O ll] 5007 A 


3.1% 


12.0% 


4.0% 


1.1% 


[S III] 9532 A 


5.8% 


5.0% 


2.0% 


2.0% 



required to achieve a given degree of convergence. The reason for 
this effect is that in a softer radiation field case the number of en- 
ergy packets emitted at wavelengths shorter than the Lyman limit 
will be less than in the case of a harder radiation field. A larger 
total number of energy packets then needs to be used in order to 
obtain a number of ionizing photons adequate to properly sample 
the nebula. The aim of Table|2|is merely to provide some general 
guidelines for selecting the appropriate number of energy packets 
for a particular simulation; however, as stated before, the optimum 
number should be determined for each given model, particularly in 
non-spherically symmetric cases. 

3.2 Benchmark Results 

The Lexington/Meudon Standard H II region model (HII40) was 
the first benchmark to be run and some very preliminary results 
have alr eady been presented, at the November 2000 Lexington 
meeting lErcolanol l200ll IPeauignot et alJ l20011. However, those 
results were produced when Mocassin was still under development 
and should therefore only be considered as a snapshot of the code 
at that particular stage. The code has evolved considerably since 
the November 2000 Lexington meeting and the newer results are 
presented in this section (see Tables- 
Table [3] shows the results of a comparison between the line 
fluxes obtained by Mocassin using the formal solution method and 
those obtained using the Monte Carlo method (see Section l2~8l for 
some of the more significant lines in the benchmark cases. It is clear 
that the results shown agree well; however, as expected, larger dis- 
crepancies were found for the weaker lines, whose lower numbers 
of energy packets yield lower accuracy statistics. 



3.2.1 The HII40 benchmark 

Mocassin scored eight ifs > 1.01 for the HII40 benchmark model 
(Table |8j; only three of these, however, had values greater than 
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Table 4. Lexington 2000 Standard H II region (HII40) benchmark case reults. 



Line 


Median 


GF 


HN 


DP 


TK 


PH 


RS 


RR 


BE 




















3-D 


1-D 


H/i/lfl 37 pro/« 


2 05 


2 06 


2 02 


2 02 


2 10 


2 05 


2 07 


2 05 


2.02 


2.10 






1 00 


1 00 


1.00 


1 00 


1 00 


1 00 


1 00 


1.00 


1.00 


Hp T SX76 

I 1 L 1 JO/U 


1 16 


1 19 


112 


113 


1 16 


1 18 


116 




0.114 


0.112 


C nl 717S-I- 


144 


157 


141 


U. 1J7 


110 


166 


OQft 


1 78 
u. 1 / 


0.148 


0.126 


P TT 1 IIS 


0.082 


100 


078 


094 


0.004 


085 


010 




0.082 


0.084 


P ml 1 Q07-u1 QOQ 


u.u /u 


07 1 
u.u / 1 


076 


n o^Q 


091 


060 


066 


074 


0.041 


0.041 


FN TTl 1 9? //m 

is 11 1 Z.Z.. 


0.034 


027 


037 


034 




032 


035 


030 


0.036 


0.034 


\N nl 6SX4.-t-6^4.8 

1> 11 UJ OtTUJtO 


730 


0.669 


817 


725 


0.69 


736 


723 


807 


0.852 


0.786 


[IN 1 1 J J / J J 


0054 


0050 


0054 


0050 




0064 


0050 


0068 


.0061 


.0054 


TlSJ TTTl ^7 1 11m 
11±J J /.J ^tlll 


292 


306 


261 


n 3 1 1 




292 


273 


101 

V. Jul 


0.223 


0.229 


TO il 6100-4-6161 

IU 1J UJUUTUJUj 


0086 


0094 


0086 


0088 


012 


0059 


0070 




.0065 


.0080 


TO ttI 7190-1-7110 

[W 1 1 J /jXUt/JJU 


U.UX7 


n o?Q 




03 1 


023 


032 


024 


036 


0.025 


0.022 


TO ttI ^79fi+^79Q 


2 03 


1 94 


2 17 


2 12 


1 .6 


2 19 


1 88 


2 26 


1.92 


1.75 


\r\ tttI 5 1 e //rn 

w 111 Ji.o ixiii 


1.06 


1 23 


1 04 


1.03 


0.99 


1 09 


1 06 


1 08 


1.06 


1.09 


TO TTll 1 11m 
[w 111J oo.j ^J.111 


1 22 


1 12 


1 06 


1 23 


1 18 


1 25 


1 23 


1.25 


1.22 


1.26 


TO tt il S007-4-4.QSQ 

[W 111J JUU/T47J7 


2 18 


2 21 


2 38 


2 20 


3 27 


1 93 


2 17 


2 08 


1.64 


1.70 


[W 111J HJUJ 


0037 


0035 


0046 


0041 


0070 


0032 


0040 


0035 


.0022 


.0023 


TO Tvl 75 Q nm 
iv-/ 1 v j _ — i . y fJiLii 


.0010 


.0010 


.0010 


.0010 


.0013 


.0013 


.0010 




.0010 


.0010 


[Np nl 17 8 ;/m 

1>C 11 l — .O 1A111 


195 


177 


195 


192 




181 


217 


196 


0.212 


0.209 


[Ne III] 15.5 /im 


0.322 


0.294 


0.264 


0.270 


0.35 


0.429 


0.350 


0.417 


0.267 


0.269 


[Ne III] 3869+3968 


0.085 


0.084 


0.087 


0.071 


0.092 


0.087 


0.083 


0.086 


0.053 


0.055 


[S 11] 6716+6731 


0.147 


0.137 


0.166 


0.153 


0.315 


0.155 


0.133 


0.130 


0.141 


0.138 


[S 11] 4068+4076 


.0080 


.0093 


.0090 


.0100 


.026 


.0070 


.005 


.0060 


.0060 


.0057 


[S III] 18.7 fim 


0.577 


0.627 


0.750 


0.726 


0.535 


0.556 


0.567 


0.580 


0.574 


0.569 


[S III] 33.6 fim 


0.937 


1.24 


1.43 


1.36 


0.86 


0.892 


0.910 


0.936 


0.938 


0.932 


[S in] 9532+9069 


1.22 


1.13 


1.19 


1.16 


1.25 


1.23 


1.25 


1.28 


1.21 


1.19 


[S iv] 10.5 (im 


0.359 


0.176 


0.152 


0.185 


0.56 


0.416 


0.388 


0.330 


0.533 


0.539 


10 3 x A(BC 3645)/A 


5.00 


4.88 




4.95 




5.00 


5.70 




5.47 


5.45 


^Inncr / K 


7653 


7719 


7668 


7663 


8318 


7440 


7644 


7399 


7370 


7480 


<T[N P N C ] >/K 


8026 


7940 


7936 


8082 


8199 


8030 


8022 


8060 


7720 


7722 


iW10 19 cm 


1.46 


1.46 


1.46 


1.46 


1.45 


1.46 


1.47 


1.46 


1.46 


1.49 


<He+ > / <H+ > 


0.767 


0.787 


0.727 


0.754 


0.77 


0.764 


0.804 


0.829 


0.715 


0.686 



GF: G. Ferland's Cloudy; PH: J.P Harrington code; DP: D. Pequignot's Nebu; TK: T. Kallman's XStar; HN: H. Netzer's 
Ion; RS: R. Sutherland's Mappings; RR: R. Rubin's Nebula; BE: B. Ercolano's Mocassin. 



1.3. Amongst these, ifs > 1.10 are obtained for the [O III] 
5007+4959 (if = 1.18) and for [O III] 4363 (if = 1.45); the ra- 
tio of these lines is often use d as a temperature diagnostic (see, 
for example. lOsterbrockll 19891 pages 119-125). Mocassin predicts 
jA49w+jA5oo7 _ 745 4 thi s value is higher than the value obtained 
by the other codes, in fact median value obtained for the ratio of 
these line fluxes by the other codes is equal to 589.2. This is fully 
consistent with Mocassin predicting a slightly lower temperature 
(if= 1.027) for this benchmark than do the other codes. 

Finally, the number of median values obtained for this bench- 
mark case is ten, which compares very well with the other codes' 
median scores, ranging from three to ten (see Table|9j- 



3.2.2 The HII20 benchmark 

Mocassin scored seven ifs for the low excitation H II region 
(HII20) benchmark model. None of these, however, has a value 
greater than 1.3. As in the HII40 nechmark case, the mean tem- 
perature, weighted by N P N C , predicted by Mocassin for this model 
is also slightly lower (;/= 1.034) than the other models' predictions. 

Five median values were obtained by Mocassin for this bench- 
mark case, while the other codes scored between three and eleven 
(see Tabled- 



3.2.3 The PN150 benchmark 

The optically thick high-excitation planetary nebula (PN150) is 
the most demanding of the benchmark cases in terms of physi- 
cal processes and atomic data required. Mocassin's score for this 
model was very good (Table |8j, obtaining only six ifs, with none 
of those being higher than 1.3 and only one slightly higher than 
1.1 (C II A1335, if = 1.13). As has already been discussed by 
IPeauignot et alj feoOll. there seems to be a dichotomy between the 
GF, HN and DP models (and, now, also the BE model) on the one 
hand, and the TK, PH and RS models on the other. The former 
group obtained very few ifs largest than 1.1, indicative of a tighter 
agreement. This coherence can probably be attributed to a more re- 
cent updating of atomic data and a more careful treatment of the 
diffuse radiation field transfer. These four codes also obtained a 
larger H/3 flux for this model, which can probably be ascribed to 
secondary photons from heavy ions. PH is the only classical code 
here with a fully iterative spherically symmetric radiative transfer 
treatment (since RR only computed H II regions); this could also 
be the reason for the relatively larger number of if s scored by the 
PH code for this model. 

The score for median values obtained by Mocassin for the 
PN150 optically thick planetary nebula is extremely good, obtain- 
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Figure 1. Electron temperature (top panels) and the fractional ionic abundances of oxygen (middle panels) and carbon (bottom panels), as a function of nebular 
radius, for the H II region benchmark cases HII40 (left-hand panels) and HII20 (right-hand panels). 
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Figure 2. Electron tempearture (top panels) and the fractional ionic abundances of oxygen (middle panels) and carbon (bottom panels), as a function of nebular 
radius, for the planetary nebula benchmark cases PN150 (left-hand panels) and PN75 (right-hand panels). 
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Table 5. Lexington 2000 low excitation H II region (HII20) Benchmark case results. 



Line 


Med 


GF 


HN 


DP 


TK 


PH 


RS 


RR 


BE 




















3-D 


1-D 


Hfl/10 36 erp/s 


4.91 


4.85 


4.85 


4.83 


4.9 


4.93 


5.04 


4.89 


4.97 


5.09 






1 00 


1 00 


1 00 


1 00 


1 00 


1 00 


1 00 


1.00 


1.00 


Hp T S876 
nc i jo/u 


0074 


0072 


008 


0073 


008 


0074 


01 10 




.0065 


.0074 


l_ 11 J ZjZJt 


04.6 


054 


047 


046 


040 


0^0 


o o^s 


063 


0.042 


0.031 


FN ttI 1 99 //m 

1> 11 1Z.Z..U111 


071 


068 




072 


007 


072 


071 


071 


0.071 


0.070 


[NT ttI 6SR4.4-6^4.X 

1> 11 UJOt'TUJT'O 


823 


745 


786 


785 


925 


0.843 


803 


915 


0.846 


0.771 


fN ttI 


0028 


0028 


0024 


0023 


0029 


.0033 


0030 


0033 


.0025 


.0021 


fN tttI S7 % nm 

1> 111J *J ' --J Ltlll 


0030 


0040 


0030 


0032 


0047 


0031 


0020 


0022 


.0019 


.0032 


TO ll 6300+6363 


.0060 


.0080 


.0060 


.0063 


.0059 


.0047 


.0050 




.0088 


.0015 


TO nl 71904-7H0 
\\J 111 / j^Ut/ JjU 


0086 


0087 


0085 


0089 


0037 


0103 


0080 


1 00 


.0064 


.0051 


TO nl 179(S4-179Q 


1 09 


1 01 


1 13 


1 10 


1 04 


1 22 


1 08 


1 17 


0.909 


0.801 


rn tttI si r /*m 

IVJ 111J Jl.O ^illl 


001 9 


0014 




001 9 


001 f\ 


001 ^ 


00 1 


0008 


.0010 


.0011 


[0 III] 88.3 


.0014 


.0016 


.0014 


.0014 


.0024 


.0014 


.0010 


.0009 


.0012 


.0013 


[0 III] 5007+4959 


.0014 


.0021 


.0016 


.0015 


.0024 


.0014 


.0010 


.0010 


.0011 


.0012 


[Ne II] 12.8 /im 


0.273 


0.264 


0.267 


0.276 


0.27 


0.271 


0.286 


0.290 


0.295 


0.296 


[S II] 6716+6731 


0.489 


0.499 


0.473 


0.459 


1.02 


0.555 


0.435 


0.492 


0.486 


0.345 


[S ll] 4068+4076 


0.017 


0.022 


0.017 


0.020 


0.052 


0.017 


0.012 


0.015 


0.013 


.0082 


[S III] 18.7 ^tm 


0.386 


0.445 


0.460 


0.441 


0.34 


0.365 


0.398 


0.374 


0.371 


0.413 


[S in] 33.6 ^tm 


0.658 


0.912 


0.928 


0.845 


0.58 


0.601 


0.655 


0.622 


0.630 


0.702 


[S in] 9532+9069 


0.537 


0.501 


0.480 


0.465 


0.56 


0.549 


0.604 


0.551 


0.526 


.582 


10 3 x A(BC 3645)/A 


5.57 


5.54 




5.62 




5.57 


5.50 




6.18 


6.15 




6765 


7224 


6815 


6789 


6607 


6742 


6900 


6708 


6562 


6662 


<T[N P N C ] >/K 


6662 


6680 


6650 


6626 


6662 


6749 


6663 


6679 


6402 


6287 


K O ut/10 18 cm 


8.89 


8.89 


8.88 


8.88 


8.7 


8.95 


9.01 


8.92 


8.89 


8.92 


<He+ > / <H+ > 


0.048 


0.048 


0.051 


0.049 


0.048 


0.044 


0.077 


0.034 


0.041 


0.048 



GF: G. Ferland's Cloudy; PH: J.P Harrington code; DP: D. Pequignot's Nebu; TK: T. Kallman's XStar; HN: H. Netzer's 
Ion; RS: R. Sutherland's Mappings; RR: R. Rubin's Nebula; BE: B. Ercolano's Mocassin. 



ing the highest value of fifteen medians, above the other codes 
which obtained between two and thirteen (see Table|9j- 



3.2.4 The PN75 benchmark 

The optically thin planetary nebula (PN75) benchmark model is not 
a radiation bounded case, but a matter bounded one and, in fact, the 
outer radius is given as an input parameter to all codes and fixed at 
7.5 x 10 19 cm. For this reason, for this particular model there is not 
a straightforward conservation law for the absolute flux of H/3. This 
can be used to explain, in part at least, the relatively poor scores of 
the GF code for low ifs (Table[S}, since, for one reason or another, 
its predicted H/3 flux deviated somewhat from the median value, 
thus shifting all the other line intensities, given in H/9 units. The 
PH code also obtained an H/9 flux which deviated from the median 
value; in this case, however, the number of total if s stayed low (=5) 
and no if > 1.30 was obtained. Mocassin, however, obtained a low 
number of ifs for this relatively difficult case, scoring nine ifs in 
total, with none of those having a value greater than 1.30. 

Mocassin obtained a score of thirteen median values for this 
benchmark case, which compares well with the scores obtained by 
the other codes for this benchmark, ranging from eight to eighteen 
median values. 



4 DISCUSSION 

The overall performance of Mocassin for the four benchmarks was 
very satisfactory, as shown by Table [8] The results obtained from 



the one-dimensional mode of Mocassin are, in general, in very good 
agreement with those obtained using the fully three-dimensional 
Mocassin models. One noticable difference, common to all four 
benchmarks, is that the kinetic temperature at the illuminated in- 
ner edge of the nebula, Tinner, is higher for the one-dimensional 
Mocassin results and closer to the values obtained by the other one- 
dimensional codes included in the comparison. This is an obvious 
effect caused by the coarseness of the three-dimensional grid: since 
all the physical properties of the gas are constant within each vol- 
ume element, then the electron temperature of a given cell will be 
mainly representative of the kinetic temperature at its centre. From 
this, it naturally follows that the coarser the grid is, and the larger 
the cells, then the further the kinetic temperature at the centres of 
the cells adjacent to the inner radius will be from the true value at 
the inner radius. 

The electron temperatures, < T[7V P 7Y C ] > and Tinner, pre- 
dicted by Mocassin for the Lexington benchmark models tend, in 
particular in the H II region cases, towards the lower limit of the 
scatter. In the case of Tinner, this seems to be a characteristic of all 
co des using an exact treat ment for the radiative transfer. As noted 
by Pequien ot et alj 1200 ll) . the kinetic temperatures calculated by 
codes with exact transfer tend to be lower in the innermost layers 
of the nebula, as the ionizing radiation field there is softer. Only 
two codes in the Lexington benchmarks treated the radiative trans- 
fer exactly, namely Rubin's Nebula (RR) and the Harrington code 
(PH) and, in fact, Mocassin's results for the kinetic temperatures 
generally agree better with those two codes' predictions. For the 
standard H II region benchmark (HII40), Mocassin's kinetic tem- 
perature at the inner edge of the nebula, Ti nner , agrees extremely 
well with the predictions of the RR and PH codes. Similar results 
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Table 6. Lexington 2000 thick planetary nebula (PN150) benchmark case results. 



Line 


Med 


GF 


HN 


DP 


TK 


PH 


RS 


BE 

3-D 


1-D 


Ha/10 35 ercrA; 


2.79 


2.86 


2.83 


2.84 


2.47 


2.68 


2.64 


2.79 


2.89 


U., 4.R61 




1 no 

1 .uu 


1 00 


1 00 


1 00 


1 00 


1 00 


1 00 


1 00 


Hp t 5R76 

1LC 1 JO/U 


0.104 


0.110 


0.129 


118 


096 


096 


095 


104 


1 06 


Hp tt 4686 

1LC 11 lUOU 


328 


324 


304 


305 


341 


333 




333 


320 


111 ZjZjt 


293 


277 


277 


293 


346 


450 


141 


339 


330 


C II 1335 


119 


121 


116 


130 




119 




103 


104 


r tttI 1907+1909 

111] 1W/T17U7 


0.174 


1.68 


1.74 


1.86 


1.69 


1.74 


1.89 


1 72 


1 71 


C TV 1 549+ 


2.16 


2.14 


2.43 


2.16 


0.154 


2.09 


3.12 


2.71 


2.65 


TN t1 5900+5 19R 

[IN 1J JZUV^TJOO 


0.012 


0.013 


0.022 


0.010 




0.020 


0.005 


.0067 


0.012 


risj ttI 65R44-654R 
[IN 11J UJoH-tUJH-o 


1 17 


1 1 s 


1 1 6 


1 1 8 
1.10 


1 01 


1 3S 


1 1 7 
1.1/ 


1 43 


1 37 


[W ttI SISS 
[IN 1 1 J J / JJ 


017 


017 


0.016 


017 


0.020 


0.023 


0.016 


0.022 


.0021 


1ST tttI 1 74.Q4. 
IN 111 I 1 / 17T 


111 


106 


109 


132 


184 


1 3Q 

U. 1J7 


091 


111 


110 


[N tttI 57 1 itm 

IN 111 ~J 1 --J /.Mil 


129 


129 


133 


134 


12 


135 


126 


120 


122 


IN 1 V I 1 H-O / T 


O 1 68 


n 1 qq 

KJ. 1 yy 


178 


192 


O 1 S4 

KJ. 1 J+ 


141 


1 68 

KJ. IOO 


162 


O 1 SQ 
u. ljy 


N V 1 740+ 

IN V 1 /.HUT 


0.147 


0.147 


159 


154 


0.055 


107 


248 


147 


0.145 


rn ii 6i 1 //m 

[W 1J UJ.l ^JJ.11 


0.020 


0.024 


017 


0.025 




.0072 


0.049 


0.010 


0.01 1 


ro t1 6ioo+6i6i 


0.135 


0.144 


0.126 


0.135 


0.245 


0.104 


0.101 


0.163 


0.153 


[O ttI 3726+3729 


2.11 


2.03 


1.96 


2.32 


2.11 


2.66 


1.75 


2.24 


2.25 


TO TTTl 51 8 ;/m 

V J 111 Jl.O I_tlll 


1.39 


1.30 


1.45 


1.42 


0.954 


1.39 


1.28 


1.50 


1.52 


fn ml RR 1 //m 

111 OO. J Ltlll 


274 


0.261 


0.292 


0.291 


27 


0.274 


0.252 


0.296 


0.299 


TO tttI 5007+4959 
\kj 111 ~j\.jyi 1 ^^y .j y 


21 4 


21 4 


22 2 


21.1 


26 


20.8 


16.8 


22.63 


22 52 


rn tttI 4i6i 

KJ 111 lJUJ 


155 


152 


151 


156 


249 


155 


109 


169 


166 


fn TVl 75 9 /;m 


3.78 


3.45 


3.16 


3.78 


3.95 


4.20 


4.05 


3.68 


3.60 


n TVl 1401+ 
U IV 1 ItUJT 


2 30 


183 


236 


324 


0.357 


0.225 




203 


201 


n vl 121 h+ 

\.y V 1^1 G ] 


0.169 


0.165 


0.189 


0.170 


0.142 


0.097 


0.213 


0.169 


0.168 


n vt 1 014+ 

V J VI IV^JT^T^ 


0.025 


0.028 


0.026 


0.022 


0.026 


0.014 




0.025 


0.026 


[Np ttI 19 R /im 

[INC 11J 1Z.O pJ.ll 


0.030 


0.028 


0.032 


0.030 


0.020 


0.027 


0.043 


0.030 


0.031 


[TnTp tttI 1 S S / <m 
[INC 111J 1J.J fJjlll 


1.97 


1.88 


1.97 


1.92 


1.73 


2.76 


2.71 


2.02 


2.03 


[Np tttI 1R69j-196R 

[INc 111J J0O7TJ7UO 


9 61 


2 64 


2 32 


2 25 


9 86 
z.ov 


3 04 


9 S6 


9 63 

i.OJ 


2 61 


ri\Tp tvI 949^4- 

[1NC 1VJ 


O 791 


n 707 
U. ikj 1 


712 


O 78S 

KJ. t OJ 


1 1 3 
1.1 j 


O 793 

U. / L.J 


O 8^9 


O 74Q 

KJ. 1 17 


O 741 

KJ. Ill 


[Np vl 1496+1146 

INC V JtZUTJJlU 


692 


721 


0.706 


0.661 


1 07 


583 


591 


692 


0.687 


[Np vl 94 9 11m 

INC V Z.t.Z. Ltlll 


980 


997 


98 


928 


1 96 


936 


195 


1 007 


0.997 


[1\Tp VTl 7 61 /(TTI 

[INC V 1J /.Oj ^illl 


O 076 
u.u / 


n 1 07 

KJ. IKJ 1 


n 07 s 

U.U / J 


077 
u.u / / 


692 


011 




050 

U.UJU 


OS 1 

U.U J 1 


Mct TT 9798+ 

IVlii 11 Z. / 70T 


1.22 


2.22 


2.10 


1.22 


0.023 


0.555 


0.863 


2.32 


2.32 


[IVlg 1VJ H.HV f_tlll 


111 


121 


111 


107 
KJ. 1U / 


1 3 

KJ. 1J 


O 049 

KJ.KJ^A 


K). 1 1 J 


111 


10Q 

KJ. lKJy 


\AAo~ vl 5 61 / »m 
[IVlg V J J.U1 fJjlll 


0.144 


0.070 


0.132 


0.162 


0.18 


0.066 




0.156 


0.156 


re; ttI Q 

[Ol 11J Jt.O £illl 


O 1 68 


O 1 5S 

KJ. 1JJ 


1 68 

KJ. 1UO 


O 1 5Q 
u. ljy 


O 963 


O 9S3 


O 1 30 

U. 1JU 


9S0 

U.ZJU 


O 963 

KJ.ZKJ3 


Si II] 2335+ 


0.159 


0.160 


0.155 


0.158 


0.20 




0.127 


0.160 


0.164 


Si III] 1892+ 


0.382 


0.446 


0.547 


0.475 


0.321 


0.382 


0.083 


0.325 


0.316 


Si IV 1397+ 


0.172 


0.183 


0.218 


0.169 


0.015 


0.172 


0.122 


0.214 


0.207 


[S 11] 6716+6731 


0.370 


0.359 


0.37 


0.399 


0.415 


0.451 


0.322 


0.357 


0.370 


[S 11] 4069+4076 


0.077 


0.073 


0.078 


0.086 


0.19 


0.077 


0.050 


0.064 


0.063 


[S III] 18.7 


0.578 


0.713 


0.788 


0.728 


0.15 


0.488 


0.578 


0.495 


0.505 


[S III] 33.6 fim 


0.240 


0.281 


0.289 


0.268 


0.06 


0.206 


0.240 


0.210 


0.214 


[S III] 9532+9069 


1.96 


2.07 


2.07 


1.96 


0.61 


1.90 


2.04 


1.89 


1.92 


[S IV] 10.5 /jm 


2.22 


2.09 


1.65 


1.76 


2.59 


2.22 


2.25 


2.25 


2.22 


^"innerV K 


18100 


18120 


17950 


18100 


19050 


17360 


19100 


16670 


17703 


<T[N P N C ] >/K 


12110 


12080 


13410 


12060 


13420 


12110 


11890 


12150 


12108 


i? ou t/10 17 cm 


4.04 


4.04 


3.90 


4.11 


4.07 


4.04 


3.98 


4.11 


4.19 


<He+ > / <H+ > 


0.704 


0.702 


0.726 


0.714 


0.79 


0.696 


0.652 


0.702 


0.711 



GF: G. Ferland's Cloudy; PH: J.P Harrington code; DP: D. Pequignot's Nebu; TK: T. Kallman's XStar; HN: H. 
Netzer's Ion; RS: R. Sutherland's Mappings; BE: B. Ercolano's Mocassin. 



are obtained for the low excitation H II region benchmark (HII20), 
where, again, Mocassin's Ti nner agrees with the results of PH and 
RR. In both H II regions benchmark cases, however, Mocassin pre- 
dicted a value which was about 250 K lower than the median for 
< T[N P N C ] >, obtaining an if = 1.027 for the HII40 case and 
if=l .034 for the HII20 case. The cause of this small discrepancy 
is not clear to us. 



Unfortunately, R. Rubin's code, Nebula, was not designed 
to treat planetary nebulae and, therefore, the only exact one- 
dimensional radiative transfer code available for the optically thick 
planetary nebula (PN150) and the optically thin planetary nebula 
(PN75) benchmarks is the Harrington code (PH). For PN150, Mo- 
cassin's Ti„„er is in reasonable agreement with PH's prediction, 
particularly if the prediction from the one-dimensional Mocassin 
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Table 7. Lexington 2000 optically thin planetary nebula (PN75) benchmark case results. 



Line 


Med 


GF 


HN 


DP 


PH 


RS 


BE 
















3-D 


1-D 




5.71 


6.08 


5.56 


5.74 


5.96 


5.69 


5.65 


5.63 


Wa 486 1 




1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


JT1C 1 JO/U 


131 


130 


144 


132 


126 


0.125 


0.132 


0.132 


Hp TT 4686 

JTT.C 11 WOU 


087 


085 


089 


087 


087 




0.093 


0.094 


V 111 Z.JZ.JT 


0.039 


0.023 


0.047 


0.040 


0.044 


0.034 


0.038 


0.043 


r n 1335 

V 11 UJ.J 


0.089 


0.096 


0.089 


0.101 


0.085 




0.086 


0.085 


r ml 1907+1909 

v. nil iy\j i ~ Lyy/y 


0.790 


0.584 


0.96 


0.882 


0.602 


1.00 


0.698 


0.709 


C TV 1 549+ 

IV 1JT7T 


0.354 


0.298 


0.480 


0.393 


0.291 


0.315 


0.414 


0.463 


IN ill 6584+6548 

[IN 11J UJOtTTjJtO 


0.098 


0.069 


0.097 


0.089 


0.108 


0.119 


0.100 


0.087 


N Til 5755 

IN 11 J / JJ 


.0012 




.001 1 


.0012 


0013 


.0020 


.0011 


.0010 


N tttI 174Q+ 

1 N 111 1 / i7T 


043 


029 


059 


056 


038 


0.048 


0.038 


0.039 


[N tttI S7 ^ //m 

IN 111 ~J 1 .-J U/lll 


397 


371 


0.405 


0.404 


390 


405 


0.336 


0.334 


N ivl 1487+ 


0.018 


0.019 


0.024 


0.020 


0.012 


0.01 1 


0.017 


0.020 


TO ttI ^7?64-^7?Q 

IVJ 1 1 J D 1 ZUTJ / 


262 


178 


262 


266 


262 


311 


0.234 


0.205 


[O tttI 5007+49^9 


11.35 


10.1 


13.2 


11.7 


10.1 


11.8 


11.0 


11.1 


[n ml 43<S3 

W H1J " 'J 


0.060 


0.046 


0.077 
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0.056 
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rn tttI S1 8 //m 


1.98 


1.94 


2.09 


1.94 


1.95 


2.02 


2.07 


2.07 
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V 7 111 OO.J 1_£111 


1.12 
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1.13 


1.12 


1.07 


1.12 


1.14 


1.14 


rn tvI 9s q nm 

\KJ ivl ZJ.7 fJilll 


814 


n 7^7 


741 


859 


821 


n 807 


0.894 


0.942 


O TVl 14(Y*4- 

V 7 1 V l^UJT 


013 


009 


015 


014 


093 




0.013 


0.015 


TlSTp Til 19 8 nm 

INC 11 1Z..O U/lll 


012 


012 


012 


012 


012 


017 


0.013 


0.012 


[KTp tttI 1 S S / <m 
[INC 1 1 1J 1 J.J fJjlll 


948 


883 


95 


902 


1 32 


1 35 


0.946 


0.949 


\Ne TTTl 1869+1968 

[INC 1 1 1J JOU7TJ7UO 


872 


784 


948 


818 


919 


1 10 


0.826 


0.838 


[Ne TVl 2421+ 


0.030 


0.028 


0.032 


0.036 


0.027 


0.020 


0.034 


0.039 


Ma TT 7798+ 


0.102 


0.086 


0.14 


0.1 1 1 


0.07 1 


0.093 


0.114 


0.106 


[A/Tct tvI 4 49 ;/m 


.0062 


.0021 


.006 


.0075 


.0065 


.0050 


.0068 


.0072 


[Si TTl 14 8 //m 
yJi ill J".u Ltiii 


0.029 


0.025 


0.034 


0.025 


0.060 


0.004 


0.061 


0.052 


Cj ttI 7IIS4. 

Ol 11J iJJJT 


.0057 


.0037 


.0078 


.0054 




.0010 


.0062 


.0052 


Si III] 1892+ 


0.104 


0.087 


0.16 


0.136 


0.101 


0.019 


0.107 


0.110 


Si IV 1397+ 


0.017 


0.017 


0.023 


0.018 


0.013 


0.023 


0.016 


0.018 


[S II] 6716+6731 


.0020 


0.023 


0.036 


0.029 


0.013 


0.016 


0.017 


0.013 


[S II] 4069+4076 


.0017 


.0022 


.0034 


.0030 


.0013 


.0010 


.0012 


.0010 


[S III] 18.7 


0.486 


0.619 


0.715 


0.631 


0.316 


0.357 


0.285 


0.266 


[S in] 33.6 


0.533 


0.702 


0.768 


0.684 


0.339 


0.383 


0.306 


0.285 


[S III] 9532+9069 


1.20 


1.31 


1.51 


1.33 


0.915 


1.09 


0.831 


0.777 


[S IV] 10.5 /nm 


1.94 


1.71 


1.57 


1.72 


2.17 


2.33 


2.79 


2.87 


10 3 x A(BC 3645)/A 


4.35 


4.25 




4.25 


4.35 


4.90 


4.54 


4.56 


^inner/ K 


14300 


14450 


14640 


14680 


14150 


13620 


14100 


14990 


<T[NpN e ] >/K 


10425 


9885 


11260 


10510 


10340 


10510 


10220 


10263 


iJout/lO^cm 




7.50 


7.50 


7.50 


7.50 


7.50 


7.50 


7.50 


<He+ > / <H+ > 


0.913 


0.912 


0.92 


0.914 


0.920 


0.913 


0.911 


0.908 


r(lRyd) 


1.47 


1.35 


1.64 


1.61 


1.47 




1.15 


1.29 



GF: G. Ferland's Cloudy; PH: J.P Harrington code; DP: D. Pequignot's Nebu; HN: H. Netzer's Ion; 
RS: R. Sutherland's Mappings; BE: B. Ercolano's Mocassin. 



run is considered, since, as discussed earlier, this represents a mea- 
surement of the temperature taken closer to the inner edge of the 
nebula. The Mocassin result for < TfiVpiVe] > is within the scat- 
ter and, in particular, BE and PH agree very well for this observ- 
able. Note that only HN and TK obtain higher temperatures for this 
model; moreover, the TK computation was carried out with a new 
code, still under development, primarily designed for X-ray stud- 
ies. That code could not treat the diffuse radiation field, leading to 
problems for the hard radiation field cases, such as PN150. Finally, 
for the PN75 benchmark planetary nebula, Mocassin's T inner is 
within the scatter (the prediction from the one-dimensional model 
is actually at the higher limit of it) and in reasonable agreement 
with PH's prediction; the result for < T[N-pN e ] > is also within 
the scatter and is in very good agreement with the prediction of 



the PH code. Once again, only HN predicts a higher value for this 
quantity, while TK's results for this model are not available. 

The models presented in this chapter were all run using a 
13x13x13 grid and, since they are all spherically symmetric, the 
ionizing source was placed in a corner of the grid. The number of 
energy packets used to sample the grids and bring them to conver- 
gence varied from three to five million. As has already been dis- 
cussed, the accuracy of the results depends both on the spatial sam- 
pling (i.e. the number of grid cells) and on the number of energy 
packets used. It is clear, however, that the latter also depends on the 
number of points to be sampled, so if, for example, in a given sim- 
ulation the number of grid cells is increased from n x xn y xn z to 
n' x xn' y xn' z , then the number of energy packets used must also be 
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Table 8. Summary of isolation factors, i/'s, for the benchmark cases 
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Table 9. Summary of median values for the benchmark cases 



Case Total GF HN DP TK PH RS RR BE 



HII40 31 8 8 10 3 9 9 5 10 
HII20 24 4 7 7 3116 4 5 



SubtotHII 12 15 17 6 20 15 9 15 



PN150 49 9 11 12 2 13 4 - 15 
PN75 40 10 8 19 - 16 13 - 13 



SubtotPN 20 19 31 (2) 29 17 - 28 



increased from iVpackcts to iVp ackcts = * ^ * ^ • iVpackcts- 
However for these relatively simple cases the three-dimensional 
grid specified above was found to be sufficient to produce accept- 
able results. In fact, since the benchmark models are spherically 
symmetric then, although the number of sampling points along 
each orthogonal axis is only 13, this is the equivalent of a one- 
dimensional code with 273 radial points, which is the number of 
different values of r given by all the (x,y,z) combinations. This is 
clearly demonstated in figures \\\ to [2] where the number of data 
points and the spacing between them shows that the spatial sam- 
pling is indeed appropriate. The plots also show that the number 
of energy packets used in the simulations was sufficient, since the 



scatter of the ordinate values for a given r, which is essentially a 
measure of the error bars, is very small. The largest scatter was 
obtained in the plots for the HII20 benchmark (Figure^; this is a 
very soft ionizing radiation field case and a larger number of energy 
packets is probably required in order to reduce the scatter shown 
and increase the accuracy of the results. For the purpose of this 
benchmark exercise, however, the accuracy achieved for HII20 is 
sufficient to produce satisfactory results. 



5 CONCLUSIONS 

A fully three-dimensional photoionization code, Mocassin, has 
been developed for the modelling of photoionised nebulae, using 
Monte Carlo techniques. The stellar and diffuse radiation fields are 
treated self-consistently; moreover, Mocassin is completely inde- 
pendent of the assumed nebular geometry and is therefore ideal for 
the study of aspherical and/or inhomogeneous nebulae, or nebulae 
having one or more exciting stars at non-central locations. 

The code has been successfully benchmarked against estab- 
lished one-dimensional photoinizat ion codes for sta ndard spheri- 
cally symmetric model nebu lae (see Pecraianot 1986; Ferla nd et alJ 
ll995tlPeauignot et all200ll) . 

Mocassin is now ready for the application to real astronomical 
nebulae and it should provide an important tool for the construction 
of re alistic nebular models. A companion paper ( Ercola no et alJ 
120021 Paper II) will present detailed results from the modelling 
of the non-spherically symmetric PN NGC 3918. Resources per- 
mitting, it is intended to make the Mocassin source code publicly 
available in the near future. 
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Appendix A: Atomic Data References 

Free-b ound emission for hydrogenic ions (H I and He H): lFerlandl 
(1980) 

Free-bound emission for He I: lBrown & Matthews! i 197d) 
Two-photon emission for hy drogenic ions (H I and He II): 

iNussbaumer & SclrniutzN 1984ft 

Two-photon emission for He I: lDrake et alJll969l) 

Free-free emission for interaction between ions of nucleus Z and 

electrons : lAlled J 19731) 

Effective recombination coefficient to H I 2 2 S:|PengglbJ_^^£4|)_^ 
Effec tive recombination coefficient to He I 2 1 S: lAlmog & NetzeJ 
(1989) 

H I a nd He II recombination line emissivities: IStorev & Hummerl 
(1995) 

He I recombination line emissivities: lBenia mine^^^999|) 
Collision transition rates for H I 2 2 S - 2 2 P: losterbrockl <1989l 

page94) 

Cooling d ue to free-fr ee radiation from hydrogenic ions (H I and 
He II): bummed! 19941) 

Cooli ng due to free-free radiation from He I: iHummer & Storevl 
(1998) 

Cooling due to recombination of hydrogenic ions (H I and He II): 
lHummeJll994l) 
Cooling due to He 1 
Collision Ionization < 

Charge exchange with hydrogen: iKingdon & Fer land ( 1996 
Fits to calculate rates of radiat ive recombination for H- 
like, He-like, Li-like, Ne-like ions: IVerner et alJ {l996). Other 



I recombinatiorrjHjmTmeT_&StoreyN 1 9981) 
>n of hydrogen: brake & UlrichHl98(J) 
vvith hydrogen: Kingdon & Ferland ( 1996) 



ions of C, N, O, Ne: Iplauignot et alJ <1991l) . Fe XVII-XXIII: 
lArnaud & Raymond Jl992l). Other ions of Mg, Si, S, Ar, 
Ca, Fe, Ni: IShull & van Steenberd <1982ft . Other ions of Na, 
Al: lLandini &^^ons!gno^^ollsll ~ Ir99f3) . Other ions of F, P, 
CI, K, Ti, Cr, Mn, Co (exclu ding Ti I-II and Cr I-IV): 

lLandini & Monsignori FossiNl99lft 

Dielectronic reco mbination coefficients: INussbaumer & Storevl 
Il983lll986lll987l) 

Non relativist ic free-free Gaunt factor for hydrogenic ions: 
lHummeJll988T) 
Fits to opacity project da ta for the photoionization cross-sections 
(outer shell) : IVerner et alJll996t) 

Collision strengths, and transition probabilities to calculate colli- 
sionally excited line strengths from ions: 



C I Col lision strengths fromlPequignot & AldrovandiH 19761) : 5 S- 
3 P from iThomas & Nesbi tl ll975ft . Transition probabilities from 
INussbaumer & RuscdJl979ft . 

C II Collision strength s fromlHav es & Nussbaumei] ll984ft . Tran- 
sition probabilities from lNus7baume^&^torevI^^8ll) . 

C III collision st reng ths and tran s ition probabilities from 
iKeenan et ai]jl992ft and lFleming et al]{l996ft. 

C IV collision s trengths from |G au & Henrvl (l977). Transition 
probabilities from lWiese et alJ fl96^T 

Mg I Collision Strengths from Saraph (1986) JA JOM calcula- 
tions. Transition probabilities from Mendozi] i 19831) . 



1994 . Transi- 



Mg II c ollision strengths and transition probabilities from 
lMendozalll983l) . 

Mg IV Collision stren gths from |Birtle^&Zeirjpenl ll994ft . Tran- 
sition probabilities from Mendoza & Zeippen ( 1983). 

Mg V Collision stre ngths fromlButler 6^eipper1il994ft . 
tion probabilities from Me ndozal (1983). 

Mg VI Collision stre ngths fromlBhatia & Masoi] ll980ft . Transi- 
tion probabilities f rom lEidelsberg et^lT^98 l |) . 

Mg VII Collisio n stren gths from lAggarwall ll984aft 
and lAggarwall Jl984bft Transition probabilities from 
INussbaumer & RuscdJl979ft . 

Ne II Collision strength from [ Ba ves et al l ^985). Transition 
probabilities from lMendoza & Z^n3perTn983l). 

Ne III Collision stre ngths from Butler & Zeipper] ll994l) . Transi- 
tion probabilities from Mendoza & Zeipper ^^83l) . 

Ne IV Collision strengths from Giles ( 1981). Transition probabil- 
ities from lZeipped ll982l) . 

Ne V Collision stren gths from | Le nnon & Burke) il99ll) . Transi- 
tion probabilities from lNussbaumer ^n^uscd*^^3) . 

Ne VI Collision strengths from Butler & Sto rey (unpublished). 
Transition probabilities fromlWjese^UdJ^^W ). 



N I Collision strengths from Berrington et al 



1981). Transition 



N II Collision strengths from [Stafford et alJ 


J 1994! 


). Transition 


probabilities from Nussbaumer & Rusca ( 1979) 
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